function [results] = DFM_function3_Extended(settings,yraw,dates,vintageDate)
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here


%% Set up otions

v2struct(settings);
model = num2str(trend+SV_factors+(trend*(1-SV_factors))*(-0.5));
  

%% Main Code

%% Some Data Adjustments

% Load Data
disp('Loading Data...');

% Remove means
y=nan(size(yraw,1),size(yraw,2));


disp('Removing means from data')


means = zeros(1,size(y,2));
stdevs = ones(1,size(y,2));

for i = 1:size(y,2)

    if size(yraw,2) > 5 && i == 2 % For consumption, impose same mean as GDP
    means(1,i) = nanmean(yraw(:,1),1);
    y(:,i) = yraw(:,i)-(remove_means)*nanmean(yraw(:,1),1);        
    else
    means(1,i) = nanmean(yraw(:,i),1);
    y(:,i) = yraw(:,i)-(remove_means)*nanmean(yraw(:,i),1);
    end
    
end

% Remove Outliers
if remove_outliers == 1
    for i = 1:size(y,2)
        y(abs(y(:,i)-localmean(y(:,i),5))>4*nanstd(y(:,i)-localmean(y(:,i),5)),i)=nan;

    end
end


%% ONLY FOR EU COUNTRIES (soft6)
for j = [soft6]
y(:,j) = y(:,j)./nanstd(y(:,j));
end
%%

disp('OK');

%% Set dimensions of the State Space system

[T,n,nftot,k,~,p_final,VarTypeInd,WhereNonMonthly,nnonmonthly,nmonthly] = SetDimensionsOfStateSpace_Extended(y,trend,[],...
    WhereQuarterly,soft6);

%% Initial Conditions for the Parameters

    [phi, delta, Sigma_trend, Sigma_epsilon, rho, Sigma_v, lambda, Q_alpha_eps, Q_sigma_eps, Q_sigma_v] = InitialConditions_Extended(T,n,p_final,k,trend,[],prior_trend);
    nsave = nreps-nburn;

%% Allocate Space for Saving Draws

Draws_phi = nan(size(phi,1),size(phi,2),nsave);
Draws_delta = nan(size(delta,1),size(delta,2),nsave);
Draws_Sigma_trend = nan(trend,trend,nreps-nburn);

Draws_Sigma_epsilon = nan(1,1,T-2,nreps-nburn);
Draws_Q_sigma_eps = nan(1,1,T-2);   
    

Draws_lambda = nan(size(lambda,1),size(lambda,2),nsave);    

Draws_Sigma_v = nan(n,T-2,nsave);
Draws_Q_sigma_v = nan(T-2,n);
    
Draws_rho = nan(size(rho,1),size(rho,2),nsave);

Draws_Factors = nan(T-2,(1+trend)*(5+1),nsave);


factors_save = [];


%% Gibss Sampler

tStart = tic;

stream = RandStream.getGlobalStream;

savedState = stream.State;

for nrep = 1:nreps
    
    
     disp(nrep);
     disp(strcat('Model:',model));
     disp(datestr(vintageDate));
     tic
  

      [F,H,Hraw,Q,R,ystar,Tstar,x0,~,z,D,A,a] = MatricesKalman3_Juan_Extended(phi,delta,Sigma_epsilon,Sigma_trend,lambda,...
            rho,Sigma_v,y,trend,[],n,nftot,k,nnonmonthly,nmonthly,p_final,VarTypeInd);
        
        Sig0 = eye(k);
        
     [xtt_sample] = DurbinKoopman(ystar,z,a,x0,Sig0,D,F,H,A,Q,R,[]);

        for jj = 1:trend;

        trend_sample = xtt_sample(1:Tstar,jj);
        d_trend_sample = trend_sample-mlag2(trend_sample,1);
        d_trend_sample = d_trend_sample(5:end,:);                   % I supress a couple of observations to reduce dependency from initial conditions

        Sigma_trend(jj,jj) = DrawIW(d_trend_sample,prior_trend,1);

        end

       init_obs = nan(p,1);
       for i = 1:p
           init_obs(end-i+1,:) = xtt_sample(1,nftot*i+1+trend:nftot*i+1+trend);   
       end

           factors = [init_obs; xtt_sample(1:Tstar,trend+1:trend+1)];
           
           factors_save = [factors_save factors];
            
           [ALPHA, SIGMA, Q_alpha_eps, Q_sigma_eps, logsigma_factor, alpha] = SV_VAR_Simple(factors,var_priors,Sigma_epsilon,Q_alpha_eps,Q_sigma_eps,prior_SV,0,0);


         delta(1:1,1:size(a,2)) = ALPHA(1:size(a,2),1:1)';
         
         phi(:,1:size(ALPHA,1)-size(a,2)) = ALPHA(size(a,2)+1:end,:)';
         
         Sigma_epsilon = SIGMA;
         
         
       observables = y;
       
       select_factors = repmat([ones(1,trend) ones(1,1) zeros(1,nnonmonthly)],1,5+1);
       observables(:,WhereNonMonthly) = [nan(q,nnonmonthly); xtt_sample(:,logical(select_factors))*lambda(WhereNonMonthly,:)' + xtt_sample(:,1+trend+1:1+trend+nnonmonthly)];         
       
       
       [lambda, ~, ~] = ...
                            BaiWangEstimationSC_SV_new_Extended(schemenumber,observables,xtt_sample(:,1:1+trend), ...
                              rho,soft6,Sigma_v,lambda_priors,trend);
   
        idiosyncratic = y(q+1:end,:) - (xtt_sample(:,:)*Hraw');
        idiosyncratic(:,WhereNonMonthly) = xtt_sample(:,1+trend+1:1+trend+nnonmonthly);

                   
        [rho, Sigma_v] = drawSerialCorrelations_new_Simple(rho,idiosyncratic,Sigma_v,rho_priors);
  
        
        [Sigma_v,Q_sigma_v] = drawStochasticVolatilities_new_Simple(idiosyncratic,Sigma_v,rho,Q_sigma_v,prior_SV_sm);
            
        
        % Plots and Saving
        
        if nrep > nburn

         Draws_Factors(:,:,nrep-nburn) = xtt_sample(:,logical(select_factors));
         
         Draws_phi(:,:,nrep-nburn) = phi;
         
         Draws_Sigma_epsilon(:,:,:,nrep-nburn) = Sigma_epsilon;
         Draws_Q_sigma_eps(:,:,nrep-nburn) = Q_sigma_eps;
         
         Draws_Sigma_trend(:,:,nrep-nburn) = Sigma_trend;
                         
         Draws_lambda(:,:,nrep-nburn) = lambda;       
         Draws_rho(:,:,nrep-nburn) = rho;
         
         Draws_Q_sigma_v(nrep-nburn,:) = diag(Q_sigma_v)';
         
         Draws_Sigma_v(:,:,nrep-nburn) = Sigma_v;
         
        end
      
    
end

% Recompute using mean values

phi = nanmedian(Draws_phi(:,:,:),3);
delta =  nanmedian(Draws_delta(:,:,:),3);
lambda = nanmean(Draws_lambda(:,:,:),3);
rho = nanmedian(Draws_rho(:,:,:),3);
Sigma_trend = nanmedian(Draws_Sigma_trend(:,:,:),3);

Sigma_v =  nanmedian(Draws_Sigma_v(:,:,:),3);
Sigma_epsilon = nanmedian(Draws_Sigma_epsilon(:,:,:,:),4);


[F,H,~,Q,R,ystar,Tstar,x0,Sig0,z,D,A,a] = MatricesKalman3_Juan_Extended(phi,delta,Sigma_epsilon,Sigma_trend,lambda,...
            rho,Sigma_v,y,trend,[],n,nftot,k,nnonmonthly,nmonthly,p_final,VarTypeInd);
        x0(1) = 3;

     
[xtt, Sigtt, ~, Sigttminus1,~] = KalmanFilter(ystar,z,a,x0,Sig0,D,F,H,A,Q,R,[]); 
[xtt_s, Sig_s] = KalmanSmoother(xtt,Sigtt,Sigttminus1,F,D,a);


       annualise = 1;
       trend_percentiles = means(1) + (1+annualise*2)*(1+standardise*(stdevs(1)-1))*prctile(squeeze(Draws_Factors(:,1,:)),[5,16,50,84,95],2);


%% Calculation of the filtering uncertainty

trend_filtering_uncertainty_filtered = nan(Tstar,nsave);
trend_filtering_uncertainty_smoothed = nan(Tstar,nsave);
trend_parameter_uncertainty_filtered = nan(Tstar,nsave);
trend_parameter_uncertainty_smoothed = nan(Tstar,nsave);

cycle_filtering_uncertainty_filtered = nan(Tstar,nsave);
cycle_filtering_uncertainty_smoothed = nan(Tstar,nsave);
cycle_parameter_uncertainty_filtered = nan(Tstar,nsave);
cycle_parameter_uncertainty_smoothed = nan(Tstar,nsave);

for jj = 1:nsave
    
    disp(jj);

    phi = Draws_phi(:,:,jj);
    lambda = Draws_lambda(:,:,jj);
    rho = Draws_rho(:,:,jj);
    Sigma_epsilon = Draws_Sigma_epsilon(:,:,:,jj);
    Sigma_trend = Draws_Sigma_trend(:,:,jj);
    Sigma_v =  Draws_Sigma_v(:,:,jj);
    
    [F,H,~,Q,R,ystar,~,x0,Sig0,z,D,A,a] = MatricesKalman3_Juan_Extended(phi,delta,Sigma_epsilon,Sigma_trend,lambda,...
            rho,Sigma_v,y,trend,[],n,nftot,k,nnonmonthly,nmonthly,p_final,VarTypeInd);

     
        [xtt_jj, Sigtt_jj, ~, Sigttminus1_jj,~] = KalmanFilter(ystar,z,a,x0,Sig0,D,F,H,A,Q,R,[]); 
        [xtt_s_jj, ~] = KalmanSmoother(xtt_jj,Sigtt_jj,Sigttminus1_jj,F,D,a);


        % Trend
        
        trend_filtering_uncertainty_filtered(:,jj) = (squeeze(Sigtt_jj(1,1,:)));
        trend_filtering_uncertainty_smoothed(:,jj) = (squeeze(Sig_s(1,1,:)));
        
        trend_parameter_uncertainty_filtered(:,jj) = ((xtt_jj(:,1)-xtt(:,1)).^2);
        trend_parameter_uncertainty_smoothed(:,jj) = ((xtt_s_jj(:,1)-xtt_s(:,1)).^2);
        
        
        % Cycle
        cycle_filtering_uncertainty_filtered(:,jj) = (squeeze(Sigtt_jj(2,2,:)));
        cycle_filtering_uncertainty_smoothed(:,jj) = (squeeze(Sig_s(2,2,:)));
        
        cycle_parameter_uncertainty_filtered(:,jj) = ((xtt_jj(:,2)-xtt(:,2)).^2);
        cycle_parameter_uncertainty_smoothed(:,jj) = ((xtt_s_jj(:,2)-xtt_s(:,2)).^2);
  

end

thing = [3*sqrt(mean(trend_filtering_uncertainty_filtered,2)...
            +mean(trend_parameter_uncertainty_filtered,2))...
         3*sqrt(mean(trend_filtering_uncertainty_filtered,2)) ...
         3*sqrt(mean(trend_filtering_uncertainty_smoothed,2)...
            +mean(trend_parameter_uncertainty_smoothed,2)) ...
         3*sqrt(mean(trend_filtering_uncertainty_smoothed,2))]; ...
         
thing2 = [3*sqrt(mean(cycle_filtering_uncertainty_filtered,2)...
            +mean(cycle_parameter_uncertainty_filtered,2))...
         3*sqrt(mean(cycle_filtering_uncertainty_filtered,2)) ...
         3*sqrt(mean(cycle_filtering_uncertainty_smoothed,2)...
            +mean(cycle_parameter_uncertainty_smoothed,2)) ...
         3*sqrt(mean(cycle_filtering_uncertainty_smoothed,2))];
         
thing3 =          (trend_percentiles(:,4) - trend_percentiles(:,2))./2;
     
%% save and plot

results.Draws_lambda = Draws_lambda;
results.Draws_Factors = Draws_Factors;
results.Draws_Sigma_trend = Draws_Sigma_trend;


variance_factor = ((1-nanmean(Draws_phi(1,2,:),3))/(1+nanmean(Draws_phi(1,2,:),3)))/((1-nanmean(Draws_phi(1,2,:),3))^2-nanmean(Draws_phi(1,1,:),3)^2);
results.sd_percentiles = 3*sqrt(variance_factor*(prctile(squeeze(Draws_Sigma_epsilon(1,1,:,:)),[5,16,50,84,95],2)));


results.means = means;
results.stdevs = stdevs;
results.dates = dates;
results.xtt = xtt;
results.trend_uncertainty = thing;
results.cycle_uncertainty = thing2;
results.trend_uncertainty2 = thing3;
results.trend_percentiles = trend_percentiles;

end

